NAR Genomics and Bioinformatics

Abstract Ribosome profiling is widely used to map translation across the genome, yet current sORF discovery pipelines are constrained by their reliance on annotated protein coding regions and by treating biological replicates primarily as tools for quantification rather than as strict filters for reproducibility. This results in two major issues, the systematic exclusion of pseudogene derived loci from analysis and the generation of extremely large candidate sets in which true translation events are difficult to distinguish from noise. To overcome these limitations, a five-stage framework was developed that expands the search space to include pseudogene loci while enforcing dual biological replicate concordance as a hard binary filter. When applied to four human cell lines with matched replicates and benchmarked against existing tools, the framework produces dramatically smaller and more reliable candidate sets, reducing thousands to millions of predicted sORFs down to a small number of high confidence concordant loci. Across HeLa cell cycle conditions, a limited set of unique concordant sORFs was identified spanning multiple transcript classes including upstream ORFs, retained introns, and nonsense mediated decay transcripts, along with a pseudogene derived candidate. Notably, the pseudogene encoded sORF RPL26P19 was detected consistently in HeLa cells but missed or inconsistently identified by existing tools, while a separate pseudogene sORF from UBE3AP2 was observed specifically in normal primary cell types and absent in HeLa. Both candidates were later found to be supported by entries in OpenProt, even though this was not known during the analysis, which gives added confidence that the results are genuine and not biased by prior expectations. Overall, these findings suggest that pseudogenes are not just inactive remnants in the genome but can act as real, functional sources of translated sORFs in a cell type specific manner. At the same time, it becomes clear that most current tools are not designed to detect these signals effectively. This highlights the need for approaches that both broaden the regions being searched and apply stricter reproducibility checks, so that truly meaningful translation events can be identified more reliably.

Background Multi-omics integration has emerged as a powerful strategy to decode the molecular complexity of biological systems. However, the diversity of available methods, each designed with distinct assumptions, objectives, and computational requirements, makes method selection, usage and interpretation challenging for nonexpert users. Here we present BiomiX 3.0, an updated version of the BiomiX platform that extends its integration capabilities with four additional methods: Similarity Network Fusion (SNF), NEighborhood-based Multi-Omics clustering (NEMO), Data Integration Analysis for Biomarker discovery using Latent variable approaches for Omics studies (DIABLO), and PRAMIGO (Phenotyping netwoRk Application for Multi-omics InteGratiOn), a novel supervised heterogeneous graph transformer (HGT) introduced in this work. Results We benchmarked all five methods, MOFA, DIABLO, SNF, NEMO, and PRAMIGO, on two independent multiomics datasets derived from a Chronic Lymphocytic Leukemia (CLL) cohort comparing IGHV-mutated and unmutated patients, and a pulmonary tuberculosis (PTB) cohort versus healthy controls. Supervised methods (DIABLO, PRAMIGO) consistently achieved higher condition-specific discrimination as measured by the Adjusted Rand Index (ARI) and the Adjusted Mutual Information (AMI). In contrast, unsupervised methods (SNF, NEMO) revealed alternative patient stratifications driven by independent sources of biological variance while MOFA performed in a semi-supervised way occupies an intermediate position, capturing latent factors that explain both disease-associated and orthogonal sources of variance. Systematic gene-centric analysis of the top-ranked features prioritized by each method was supported by manual biological annotation of shared and method-specific signals. Across cohorts, we annotated 154 shared features (116 genes, 38 metabolites) and 60 method-unique features per cohort, demonstrating that no single integration strategy captures the full landscape of biologically relevant signals. In the CLL cohort, shared features spanned B-cell receptor biology, innate immune signaling, RAS/MAPK activation, and epigenetic regulation, while methodunique features revealed supervised-method-specific insights into vesicle trafficking (DIABLO), immune checkpoints (MOFA), and ncRNA regulation (PRAMIGO). In the PTB cohort, a convergent interferon/innate immune signature dominated shared features across all methods. Still, method-unique analysis uncovered DIABLO-specific acylcarnitine metabolic reprogramming, MOFA-specific restoration of lysosomal trafficking, and PRAMIGO-specific {gamma}{delta} T-cell and immunoglobulin repertoire diversity. Across both cohorts, SNF and NEMO proved useful for detecting biological and technical sources of variation that were orthogonal to the primary condition of interest. Ultimately, PRAMIGO uniquely enables the construction of heterogeneous graphs modeling cross-modal molecular interactions, uncovering epigenetic co-regulation programs in CLL and multi-omics inflammatory modules in PTB that are difficult to identify using conventional integration approaches. Conclusions BiomiX 3.0 provides a graphical user interface (GUI) multi-method integration environment that democratizes access to state-of-the-art multi-omics analysis. By combining both unsupervised and supervised integration strategies within a unified platform and introducing graph-based learning through PRAMIGO, BiomiX 3.0 enables researchers across disciplines with complementary tools to interrogate the biological sources of variation in their data, without requiring bioinformatics expertise.

Ribosome profiling experiments capture a snapshot of translation by sequencing ribosome-protected fragments (RPFs). Alignment to a well-curated transcriptome captures most RPF signal, but can miss ribosome activity when expressed exons are absent from the reference. Genome alignment of RPFs expands the detectable sequence space, but can expose ambiguity between homologous loci. Existing workflows often couple reference choice with the handling of ambiguously mapped reads, making it difficult to isolate the effects of either decision on downstream analyses. Here, we developed RiboFlow v2 by adding genome alignment and configurable handling of multimapping reads to the existing transcriptome workflow. Using 24 human ribosome profiling libraries and their matched RNA-seq data, we found that nucleotide-level coverage at matched coding-sequence positions and translation efficiency (TE) were broadly concordant between alignment strategies under stringent alignment filters. Even under these conditions, individual genes showed substantial differences in coverage and TE estimates. These differences included coverage in alternative exons absent from selected transcripts and ambiguous alignments between protein-coding genes and processed pseudogenes. Clustering genes by read-assignment patterns identified groups associated with pseudogene prevalence and highlighted shared exonic sequence among reference transcripts as a source of alignment discrepancies. RiboFlow v2 extends the analysis beyond selected reference transcripts and enables identification of genes whose ribosome coverage and abundance estimates depend on reference composition and how ambiguously mapped reads are handled.

Rapid and accurate identification of viral sequences underpins clinical diagnostics, epidemiological surveillance, and the safety testing of biological products, yet the established alignment-based methods such as BLAST are inherently closed-set, with accuracy tied to how well a query is already represented in the reference database. Recent advances in genomic foundation models (GFMs) have enabled alignment-free approaches to sequence classification. However, their performance relative to traditional methods remains incompletely characterised. We evaluated GFMs, specifically the Nucleotide Transformer (NT) and DNABERT families, for viral sequence classification using a dataset derived from the Reference Viral Database (RVDB). Sequence embeddings were generated with mean, max, and CLS pooling, classified by nearest-neighbour search in an embedding index and compared with BLAST as a reference. Matching BLAST provided a particularly stringent benchmark here, because as a closed-set method it compares each query directly against the exact reference sequences it stores, whereas a foundation model must rely on representations shaped by broad, general-purpose pre-training. On sequences of standard length, the best model, NTv3-650M, reached an accuracy of 0.93, approaching the 0.97 achieved by BLAST. On long genomic sequences, evaluated on a much smaller test set, the difference narrowed further, with NTv3-650M reaching 0.95 against 0.98 for BLAST. Robustness analyses showed that all GFMs were markedly more sensitive to base masking than BLAST. Both methods indexed the same labelled reference sequences, and BLAST compared queries with them directly at the nucleotide level, so approaching it from frozen, general-purpose representations without any parameter update is a meaningful result. NTv3-650M emerged as the most promising candidate for further optimisation and deployment in bioinformatics applications.

Selection of an appropriate coding exon is an important upstream step in exon deletion-based loss-of-function strategies, particularly for genes with multiple transcript isoforms. However, manual comparison of exon structures across multiple transcripts can be time-consuming and may overlook transcript-specific differences in exon boundaries and coding sequence (CDS) composition. Non3nExonFinder was developed as an interactive R/Shiny-based tool for identifying shared non-triplet coding exons across user-defined human or mouse NCBI RefSeq protein-coding transcripts. Non3nExonFinder evaluates exon/CDS annotations and coding contributions for each transcript, identifies common exons with identical genomic boundaries, and prioritizes candidates based on CDS-only status and coding-length divisibility by three. Non-triplet CDS-only exons that are internal in all selected transcripts are classified as high-confidence candidates. Independent reconstruction from the original RefSeq GTF files across 80 genes (40 human and 40 mouse), comprising 6,848 transcript-specific exon instances, showed complete concordance with Non3nExonFinder outputs, with exact recovery of all candidate sets. Non3nExonFinder provides a reproducible upstream framework for transcript-aware prioritization of exon-deletion targets.

Summary: Controlled modification of sequencing data is important for reproducible benchmarking, particularly when evaluating analyses affected by read fragmentation, divergent reference genomes and ancient DNA damage. SCAR introduces controlled mutations, fragmentation and position-specific damage into FASTA and FASTQ data. All features can be used separately or in combination, and empirical fragment-length and mismatch profiles can be supplied to closely reproduce the properties of targeted datasets. Validation analyses showed the successful implementation of the requested sequence changes, and their impact on read mapping and heterozygosity analyses. Availability and Implementation: SCAR is implemented in C++ and the latest code is available at https://github.com/Madshartmann1/SCAR. The version described in this article is archived at Zenodo at https://doi.org/10.5281/zenodo.22789736 , with validation material available at https://doi.org/10.5281/zenodo.21915766.

Motivation: High-throughput sequencing technologies have driven rapid growth of biological sequence databases. Public repositories must therefore rely on automated computational heuristics to screen submitted sequences for errors and low quality. For example, SILVA SSU Ref, which exploits the conserved nature of 16S and 18S rRNA sequences, applies strict algorithmic quality controls yet still accepts sequences with up to 30% of their nucleotides deviating from any previously accepted sequence. This permissiveness creates opportunities for the admission of modified sequences, such as biologically plausible sequences generated by DNA foundation models. The vulnerability of public sequence databases to becoming polluted or poisoned with modified sequences necessitates the development of methods to detect such sequences. Results: We present the first investigation, to our knowledge, of the detectability of modified sequences. We consider simple computationally-modified 16S rRNA sequences that pass the quality control inclusion criteria of the SILVA database, which we generate via random substitutions. We present classifiers that can distinguish such modified sequences from natural 16S rRNA using conserved motifs. Our best classifier achieves over 90% sensitivity and specificity on our testing set when using a 5% artificial mutation rate. One feature used in our classifiers, gapped k-mers constructed from universally conserved nucleotides, was conserved across all three domains of life despite relying on exact matches to patterns found in E. coli, advancing our understanding of conserved grammatical structure in small subunit rRNA sequences. Availability and Implementation: Our source code is available at https://github.com/rainhaworth/16S-Mutation-Classifiers.

We have developed a method to integrate transcriptomic coexpression networks across diverse experimental conditions within a single species. Our framework expands previous pannetwork approaches into a hypergraph based pannetwork. It first identifies coexpressed gene clusters within each individual dataset and extracts frequently co-expressed gene sets across multiple datasets using a frequent itemset mining algorithm. This process yields a hypergraph where each hyperedge is assigned a frequency (universality, U). We then extracted a subnetwork comprising high U hyperedges as the core network and formed modules within it. We applied our method to 106 Escherichia coli transcriptome datasets from the GEO database. The modularity of the core network peaked at a universality cutoff of 15, which was subsequently used to define it. Approximately 70% of the resulting core modules correspond to operons, and conversely, approximately 70% of all operons are covered by these core modules. We visualized the core network via an inter-modular network and analyzed core module dataset relationships using a modularity profile matrix. Based on these analyses, we successfully visualized the dynamic reorganization of the coexpression network in response to environmental changes in bacteria.

Background: Protein-coding DNA exhibits characteristic three-base periodicity. In unaligned DNA libraries, however, this signal should cancel if fragment boundaries are uniformly distributed across codon phases. Yet, we find that in multiple metagenomic libraries, base frequencies depended systematically on distance from fragment boundaries, producing a period-3 pattern. We investigated the origin of this unexpected signal. Results: The signal was observed in modern metagenomic datasets and widely across ancient-DNA datasets. We developed a new alignment-free algorithm to infer latent triplet-phase structure from nucleotide composition. The inferred shifts were strongly nonuniform at fragment boundaries, explaining the read-coordinate periodicity. We initially suspected a protocol artifact. To test an alternative, we derived a probabilistic framework linking nucleotide-dependent boundary selection to codon-phase frequencies. The framework showed that, because nucleotide composition differs among codon positions, preferential breakage at particular nucleotides biases the phases represented at fragment boundaries. Simulations confirmed this mechanism: purine-associated fragmentation produced strong period-3 coherence in a coding-dense bacterial genome, whereas no comparable genome-wide effect appeared in the human reference genome. Conclusions: The observed read-coordinate three-base periodicity does not require an intrinsic phase bias in the source DNA. Instead, sequence-dependent fragment formation or recovery can couple boundaries to codon phase. Purine-associated post-mortem fragmentation plausibly explains the widespread phenomenon in ancient-DNA libraries, although PCR- or library-specific mechanisms may contribute to it in modern datasets. Our alignment-free algorithm provides a practical method for estimating and normalizing phase shifts before searching for positional nucleotide patterns, preventing boundary-associated bias from being mistaken for an intrinsic biological pattern.

Advances in sample extraction techniques, improvement of library preparation methods and the dramatic reduction in sequencing costs have made complex, multifactorial experimental design the de facto standard in epigenomic studies. In contrast to the plethora of tools, workflows and libraries available to analyze gene expression data, consolidated workflows for chromatin accessibility data are lacking and analysis of such data often occurs on an ad hoc basis. Here, a scalable end-to-end workflow, ATACofthesnake, is presented that simplifies analysis of multifactorial chromatin accessibility data considerably and provides a newly created framework to handle time series data reproducibly.

After decades of investigation, analyses of genome organization remain challenging. This is because genome descriptions are very information-rich and usually rendered as a linear map of genomic features corresponding to an ordered list of position intervals on a nucleic acid macromolecular backbone. This formalism, while accurate, is not fully optimized to mathematically analyze all genome structural properties, in particular the relations between genomic features, nor to integrate large amounts of -omics data, such as transcriptomes, epigenomes, proteomes, meltomes or chromatin accessibility data, which are increasingly available. We therefore developed GenDiMug, a tool that uses a network representation of the genome as a backbone, providing a natural framework to integrate spatial and ontogenetic relationships between genomic features and multi-omics data. We illustrate how this generic approach can be applied to all sorts of genomes, both to visualize genomic topology at contig and chromosomal scales at base-pair resolution, and to analyze genome positional effects using various graph metrics. Using GenDiMug, we further reveal clumping of genes with similar expression levels in eukaryotic genomes and their possible circadian regulation in a photosynthetic alga.

Background: Pangenome variations graphs (PVG) are promising structures to study diversity and reduce reference bias in genomics. However, it is not currently possible include complete annotation of the PVG. In a previous work, we introduced GrAnnoT, a tool for transferring the annotation from a reference genome to a PVG, and from a PVG to its embedded genomes. However, GrAnnoT was limited to transferring one annotation to the graph, and thus could not include all genes. Results: We implemented new options to GrAnnoT, that can now handle annotations of multiple genomes, and integrate them in the graph. We also developed GrAPEx, a pipeline that identifies, extracts and annotates non-annotated regions in the graph. We tested those two tools to attest their performances. Our benchmarks revealed that both GrAnnoT and GrAPEx are reliable. GrAnnoT is now able to integrate correctly several annotations, increasing the amount of biological information available in the graph. Additionally, it offers useful options to study gene distribution in the PVG. GrAPEx allows rapid annotation of the whole PVG, with no reduction in the annotation quality. Conclusion: Overall, the two tools presented here allow to annotate thoroughly the whole PVG, and integrate biological annotation to the pangenome. It also reduces the previous reference bias, where only one annotation (often the reference one) could be included in the graph. This constitutes a precious addition to the pangenomic toolkit.

Small RNA (sRNA) sequencing is widely used to study siRNA, miRNA and piRNA pathways and to profile sRNAs derived from viruses, transposons and host transcripts. However, many sRNA workflows rely on command-line mapping and separate software for downstream analysis and visualisation. We present sRNAMap, an HTML5/JavaScript application for mapping and analysing sRNA sequencing data. sRNAMap accepts either raw FASTA/FASTQ reads or pre-aligned BAM files. For raw reads, sRNAMap uses a k-mer-seeded alignment algorithm that runs directly in the browser using parallel Web Workers where supported. It generates strand-specific per-base coverage, 5'- and 3';-end profiles, read-length distributions and nucleotide-bias summaries, with optional library-size normalisation. It also quantifies characteristic read-overlap signatures, including 5';-to-5'; overlaps associated with ping-pong-like processing. Additional modules quantify phased sRNA production, non-templated 3'- and 5'-additions, and sequence diversity. Analyses can be saved and restored as JSON or exported as self-contained HTML reports containing figures and analysis metadata. Benchmarking against Bowtie 2 under approximately matched alignment settings using a representative virus-infected arthropod small-RNA dataset showed highly concordant mapping profiles. Per-position depth showed an R^2 of approximately 0.998, while 5';- and 3';-end counts showed R^2 values of approximately 0.97. sRNAMap is available at https://github.com/rhparry/sRNAMap under the MIT licence.

Motivation: De novo genome assembly relies on accurate path selection through complex graphs, yet unitig graph branching points remain a major ambiguity source. Multiple extensions may be equally plausible from topology alone, but sequence-level consistency of competing paths remains largely unexplored. Model: We introduce RASTRO, a reference-free framework for scoring candidate paths at unitig branching points using nucleotide sequence information. RASTRO-S employs a sample-specific n-gram model to estimate extension plausibility via bits-per-base surprisal, while RASTRO-L uses the Evo2 pretrained DNA language model for autoregressive path scoring. Both operate directly on unitig graphs without reference genomes or assembled contigs. To evaluate branch decisions independently, we develop a read-based validation strategy using boundary k-mers spanning the junction between graph context and candidate extensions. Results: Experiments on HIV-1 and SARS-CoV-2 unitig graphs from Logan show sequence-based models provide meaningful match with read-supported choices. Among 20,000 sampled branching points, unique read-supported choices were identified for 15,319 HIV-1 branches (76.6%) and 19,882 SARS-CoV-2 branches (99.4%). RASTRO-S improved match from 54.1% to 75.8% for HIV-1 and from 58.6% to 76.4% for SARS-CoV-2 as context length rose from n=5 to n=10. RASTRO-L also improved with context length, reaching 49.96% and 49.25% match. Both outperformed random selection and the longest-extension greedy baseline. Conclusion: RASTRO provides an efficient, interpretable framework for evaluating ambiguous path decisions in unitig graphs using raw-read evidence rather than contig-level alignment. Results demonstrate the potential of lightweight, sample-adaptive sequence models for improving local decision-making in genome assembly workflows.

Protein-RNA interactions play a central role in many cellular processes, such as gene regulation, protein synthesis and viral infections. Although a few thousand protein-RNA complexes have been structurally characterized, they represent only a small fraction of all interactions. This lack of experimental data limits the ability of deep-learning based approaches such as AlphaFold3 to predict structures of protein-RNA interactions. In this study, we investigate protein-RNA binding specificity by enriching experimental structures with omics data. To this end, we present a scoring approach and associated web resource quantifying the agreement between experimental structures of protein-RNA interfaces and their omics-derived binding preferences, allowing us to identify interaction motif cores. Through key structural and evolutionary features, we further highlight that these motif cores correspond to important interface regions. We then leverage the dataset of protein-RNA complexes for which structural information can be combined with binding preferences in order to benchmark the ability of AlphaFold3 to predict protein-RNA interaction specificity. To this aim, we run AlphaFold3 predictions using different RNA inputs, from the exact sequence present in the experimental structure to a non-specific sequence, including sequences embedding the consensus binding motif from in vitro experiments. We show that despite good prediction quality, the sensitivity of AlphaFold3 to the exact RNA sequence used as input indicates signs of memorization. We examine some particular complexes and uncover challenges encountered by our scoring workflow and AlphaFold3, especially in the case of alternative binding modes. Altogether, this work evidences both the promise and current limitations of deep learning approaches for protein-RNA structure prediction, and provides a resource to guide their further development.

We present a lossless compression algorithm for indexing a protein database while supporting fast taxonomic classification in the method Centrifuger. The algorithm is a new scheme of the previously proposed run-block compression algorithm to reduce the size of the Ferragina-Manzini (FM) index, and it scales better with alphabet size than the original. On the RefSeq prokaryotic and viral protein sequences, Centrifuger reduces the memory footprint by over a third compared to the method Kaiju that builds on a plain FM-index, while having comparable running time. Furthermore, the compressed FM-index is lossless and can locate matches of arbitrary length, which helps Centrifuger achieve greater accuracy than Kraken2, a k-mer-based taxonomic classification method. We leverage the computational efficiency of Centrifuger to create an index of size 182 GB for classifying the reads against the full nr database that contains about 250 billion amino acid characters. Using this index, Centrifuger reveals different SARS-CoV-2 infection states and viral transcriptome profiles across human cell types from single-cell RNA-seq data.

Functionally linking enhancers to gene promoters at high resolution often entails analysis of 3C datasets, including viewpoint-based methods such as Capture-(Hi)C. Quantitative comparison across samples and experiments, together with visualization and integration of promoter interactions with other datasets are of critical importance. While existing statistical tools identify significant interactions and quantify underlying signal, systematically comparing the shape and extent of this signal across large sets of intervals remains challenging. Moreover, experimental artefacts, such as inefficient crosslinking, can affect downstream interpretation, but are difficult to detect systematically. Here we introduce PostChicago, an R/Bioconductor package that fills this gap. Building upon the widely used interaction-calling package Chicago, PostChicago provides a complete workflow for quality control, visualization, integration and quantitative analysis of promoter interactions between multiple samples, including annotation and quantification of genomic regions based on external datasets. PostChicago extends visual quantitative comparisons beyond individual viewpoints by computing normalized metaprofiles around large sets of intervals that visualize the distribution of interaction signal, enabling comparison between narrow and broad interaction patterns. We further demonstrate how PostChicago can identify experimental crosslinking artefacts. Together, PostChicago enables integrative, shape-sensitive analysis and quality control of promoter Capture-(Hi)C datasets at a scale that would otherwise require extensive custom scripting.

Nanopore direct RNA sequencing has emerged as a powerful technology for studying RNA and modifications at single-molecule resolution. Although originally developed for polyadenylated RNAs, protocol adaptations have extended its application to non-polyadenylated RNA species, including tRNAs and rRNAs. However, the lack of open and flexible (de)multiplexing solutions remains a limitation for large-scale sequencing studies. Existing community-driven approaches enable multiplexing of polyadenylated RNAs and, more recently, tRNAs, but do not fully address the challenges of targeted sequence-specific Nanopore sequencing of non-polyadenylated RNAs. Here, we present an improved sequence-specific direct RNA sequencing protocol for ribosomal RNAs, achieving enhanced capture and longer read lengths as demonstrated for human 28S rRNA. We complement this protocol with an expanded WarpDemuX software toolbox for flexible sequence-specific (de)multiplexing. We validate the approach using human cell lines carrying knockouts of different C/D box snoRNAs involved in site-specific 2'-O-ribose methylation of 28S rRNA. Our results demonstrate efficient multiplexed detection of sequence-specific RNA molecules and associated modification changes. While established here using rRNA, we validate our approach on another panel of non-polyadenylated RNA species, providing an open and scalable framework for targeted direct RNA sequencing and RNA modification analysis

Genomic variant calling is a central computational task in disease research and diagnostics, where somatic mutations such as JAK2 V617F serve as diagnostic markers in myeloid neoplasms. VarDict, originally written in Perl, and its reimplementation VarDictJava are widely used variant callers in the bioinformatics community. The Perl implementation requires run times of up to several hours; VarDictJava reduces this to a few minutes but can consume large amounts of main memory, which limits the number of jobs that can be run in parallel on a given machine. We present vardictcpp, a functionally identical reimplementation of VarDictJava in C++17. vardictcpp is approximately four times faster in single-core execution and 9-11 times faster using eight cores, while requiring up to 17 times less main memory. Variant calls are identical to those of VarDictJava on the evaluated whole-exome samples. The reduced memory footprint allows more concurrent jobs on existing hardware and thereby increases sample throughput, which is particularly relevant for clinical sequencing laboratories with fixed compute resources. vardictcpp is implemented in C++17 and freely available at https://github.com/MHH-Bioinformatics-Hematology/vardictcpp

High throughput sequencing necessitates robust computational tools to contextualize complex transcriptomic data. In this study, we analyzed RNA Seq data from peripheral blood of bacterial pneumonia patients and healthy controls (GEO GSE161731) to investigate systemic molecular responses to infection. Differential gene expression analysis identified significantly upregulated genes in pneumonia. We used GANGO for Gene Ontology (GO) term enrichment and BioFunctional for hierarchical exploration, providing an initial overview of affected biological processes. To deepen interpretation, we leveraged CURIE (www.curieai.tech ), an AI plat-form that cross-references GO terms with scientific literature. CURIE generates narrative explanations, link-ing molecular mechanisms, experimental evidence, and publications to each biological process, moving be-yond simple enrichment. Our combined analysis revealed a coordinated upregulation of pathways involved in DNA repair and replication, transcriptional activation, RNA processing, translation, ribosome biogenesis, mitochondrial respiration, and cofactor binding. By integrating these findings within a broad, literature-supported framework, we highlighted how multiple cellular compartments and processes collaborate to mount an effective response to infection-induced stress. This study underscores the potential of combining classical statistical enrichment with AI based interpretation for a holistic understanding of host-pathogen in-teractions at the transcriptomic level.